#include "header.h"

/* -------------------------------------------------------------------------- */
void iterate(int iter,int rk_stage_count,struct node* p,
          void (*FLUX)(int,struct node*,double**,double*)){
/* -------------------------------------------------------------------------- */
int i, j, k, j1;
double divf[5];
void update(int,int,struct node*,double**,double*);
void supout(int,struct node*);
void flux_derv_mirror(int,struct node*,double**,double*);

if(time_mode == GLOBAL){
	for(i=0;i<no_grid_points;i++)
		for(j=0;j<5;j++) var_old[i][j] = p[i].var[j];
	}
QDerv(p,var_old);

for(i=0;i<no_grid_points;i++){  /* Updating for each node begins here */
	switch(p[i].point_type){
		case CONST: break;

		case OUTER_BOUND: 
			//outer(i,p);
			break;

		case INTERPOLATE:
			interpolate(i, p);
			break;

		case SUP_OUT:
			//supout(i,p); break;
		case SOLID_WALL:
		case KUTTA:
		case FIRST_LAYER:
		case CARTESIAN:
		case INTERIOR:
			(*FLUX)(i,p,var_old,divf);
			update(rk_stage_count,i,p,var_old,divf);
			break;

		default:
			printf("Point type not assigned properly for node %d\n",i);
			printf("Type is %d\n", p[i].point_type);
			exit(0);
			break;
		}
	}


}

/* -------------------------------------------------------------------------- */
void update(int rkc,int point,struct node* p,double** var_old,double* divf){
/* -------------------------------------------------------------------------- */
int i;
double rkfact,var[5],con[5],cono[5],vel2, un;
void contoprim(int,double*,struct node*);
void correct_nwbc(int,struct node*,double**);

//rkfact = 1.0/(double)(rk_stage - rkc);

dt_cfl = p[point].dt;
for(i=0;i<5;i++) var[i] = p[point].var_old[i];

vel2   = var[1]*var[1] + var[2]*var[2] + var[3]*var[3]; 

con[0] = var[0];
con[1] = var[0]*var[1];
con[2] = var[0]*var[2];
con[3] = var[0]*var[3];
con[4] = var[4]/(GAMMA-1.0) + 0.5*var[0]*vel2;

vel2   = var_old[point][1]*var_old[point][1] + 
         var_old[point][2]*var_old[point][2] + 
         var_old[point][3]*var_old[point][3];

cono[0] = var_old[point][0];
cono[1] = var_old[point][0]*var_old[point][1];
cono[2] = var_old[point][0]*var_old[point][2];
cono[3] = var_old[point][0]*var_old[point][3];
cono[4] = var_old[point][4]/(GAMMA-1.0) + 0.5*var_old[point][0]*vel2;

switch(p[point].point_type){
	case KUTTA:
        for(i=0;i<5;i++)
           con[i]   = AAA[rkc]*con[i] + BBB[rkc]*(cono[i] - dt_cfl*divf[i]);
		//con[1] = con[2] = con[3] = 0.0;
        contoprim(point,con,p);
	break;

	case OUTER_BOUND: /* Not implemented yet */
		printf("iterate: Error, outer boundary not implemented.\n");
		exit(0);
	break;

	case SOLID_WALL:
        dt_cfl = wall_cfl*dt_cfl;
        for(i=0;i<5;i++)
           con[i]   = AAA[rkc]*con[i] + BBB[rkc]*(cono[i] - dt_cfl*divf[i]);
           contoprim(point,con,p);
        /* vel2   = p[point].var[1]*p[point].var[1] + p[point].var[2]*p[point].var[2];
		p[point].var[1] = sign(p[point].var[1])*sqrt(vel2); */
            /* Enforce zero normal velocity. Requires only normal vector. Idea
             * given by Anandhanarayanan */
		un = p[point].var[1]*p[point].n[0] + p[point].var[2]*p[point].n[1] 
		     + p[point].var[3]*p[point].n[2];
		for(i=0; i<3; i++) 
                    p[point].var[i+1] = p[point].var[i+1] - un*p[point].n[i];
		switch(wallmode){
			case SWBC:   break;
			case PSEUDO: break;
			case NWBC:   correct_nwbc(point,p,var_old); break;
			default: printf("No matching wall mode.\n");
			}
	break;

	case SUP_OUT:
	case INTERIOR:
	case FIRST_LAYER:
	case CARTESIAN:
        for(i=0;i<5;i++)
           con[i]   = AAA[rkc]*con[i] + BBB[rkc]*(cono[i] - dt_cfl*divf[i]);
        contoprim(point,con,p);
	break;
	
	default:
		printf("Point type is not assigned properly for node %d\n",point);
		printf("Type is %d\n", p[point].point_type);
		exit(0);
	
	}	


}
/* -------------------------------------------------------------------------- */
void contoprim(int point,double* con,struct node* p){
/* -------------------------------------------------------------------------- */
double vel2;

if(con[0] <= 0.0){  /* Check if the density is positive */
	printf("iterate: Density is negative at node %d\n",point);
	printf("x = %f   y = %f\n",p[point].x,p[point].y);
	exit(0);
	}

p[point].var[0] = con[0];
p[point].var[1] = con[1]/con[0];
p[point].var[2] = con[2]/con[0];
p[point].var[3] = con[3]/con[0];
vel2            = p[point].var[1]*p[point].var[1] + 
                  p[point].var[2]*p[point].var[2] + 
                  p[point].var[3]*p[point].var[3];
p[point].var[4] = (con[4] - 0.5*con[0]*vel2)*(GAMMA-1.0);
if(p[point].var[4] <= 0.0){
    printf("iterate: Negative pressure at node %d, point type = %d\n",
            point,p[point].point_type);
	printf("Conn size = %d\n", p[point].total_conn);
    printf("x = %f, y = %f, z = %f\n", p[point].x, p[point].y, p[point].z);
    exit(0);
    }
}
/* -------------------------------------------------------------------------- */
